'''
同频：把你的就业率预测（日/半月）聚合到与自变量一致的频率（比如 GDP 用季度）。

变换：为避免“趋势伪相关”，推荐做 YoY（季度→滞后4期；月度→滞后12期）。规模类（GDP/一二三产）用对数同比，占比用简单同比差分；目标（就业率）也做同比差分。

再计算 Pearson / Spearman 及 p 值（必要时做 FDR 校正）。

同频对齐

就业率预测是日/半月频，代码把它按季度取均值（resample("Q").mean()）。

GDP表本来就是季度；为了保险也 resample("Q").last() 一下。

然后按日期对齐，只保留两边都有的季度。

选变换方式（默认 YoY，同比）

目的：去掉共同趋势和季节，别算出“看起来很相关其实是一起上升”的伪相关。

季度的同比就是跟 上年同季度比（滞后 4 季）：

规模类（GDP、第一/二/三产业）：用对数同比 log(X_t) - log(X_{t-4})。

占比（第一/二/三产业占比）：用简单同比差分 X_t - X_{t-4}。

目标（就业率）：也用同比差分 Y_t - Y_{t-4}。

如果你把 TRANSFORM 改成 "diff"，就是一阶差分；改成 "level" 则用原始值（不推荐，容易伪相关）。

逐个变量算相关

对每个自变量（变换后的）和就业率（变换后的）：

算 Pearson（线性相关）和 Spearman（单调相关）；

给出 p 值（显著性）；

记录样本数 n_obs（太少就不靠谱）。

因为你一次性测了很多变量，顺带做了 FDR（多重检验校正），在表里给出 pearson_q（越小越稳）。

导出结果

把所有变量的相关、p 值、q 值、样本数，写到 corr_only_results.csv，并在控制台打印一份。

为什么要“先对齐再差分/同比”？
同频：季度对季度、月度对月度，才有可比性。

去趋势/去季节：宏观量很多都一起“稳步上升”或有明显季节，直接算相关容易得到“看着很高”的伪相关。同比/差分能把“长期趋势”和“固定季节”剥掉，专注于“变化跟变化”的真实关系。
'''


# -*- coding: utf-8 -*-
# Part 1: 相关性系数（先对齐→再做同比/差分→相关）
# 默认：季度YoY（GDP_季度数据.csv + data_processed.xlsx）
import os, numpy as np, pandas as pd
from scipy.stats import pearsonr, spearmanr
from statsmodels.stats.multitest import multipletests

# ===== 配置 =====
PATH_EMP = "/Users/linshangjin/25CCM/NKU-C/dataset_csv/data_processed.xlsx"   # 就业率预测（多Sheet，date, employment_rate）
PATH_X = "/Users/linshangjin/25CCM/NKU-C/dataset_csv/GDP_季度数据.csv"       # 季度GDP与产业结构
OUT_DIR   = "/Users/linshangjin/25CCM/NKU-C/t2/三次产业/相关性"
FREQ     = "Q"                                # 'Q'=季度，若是月度源就改为 'M'
TRANSFORM = "yoy"                             # 'yoy'（推荐）或 'diff' 或 'level'
OUT_CSV  = "/Users/linshangjin/25CCM/NKU-C/t2/三次产业/相关性/mv_corr_quarterly_yoy.csv"

# ===== 工具 =====
def read_employment_xlsx(path):
    xl = pd.ExcelFile(path); lst = []
    for sh in xl.sheet_names:
        t = xl.parse(sh)
        assert {"date","employment_rate"}.issubset(t.columns), f"{sh} 需含 date, employment_rate"
        t = t[["date","employment_rate"]].copy()
        t["date"] = pd.to_datetime(t["date"])
        lst.append(t)
    df = pd.concat(lst, ignore_index=True).dropna(subset=["date"]).sort_values("date")
    df["employment_rate"] = pd.to_numeric(df["employment_rate"], errors="coerce")
    return df.set_index("date").sort_index()

def read_vars_table(path):
    # 兼容 CSV/Excel，日期列名可能是 'date' 或 'Unnamed: 0'
    if path.lower().endswith(".csv"):
        d = pd.read_csv(path, encoding="utf-8")
    else:
        d = pd.read_excel(path)
    if "date" in d.columns:
        d = d.rename(columns={"date":"date"})
    elif "Unnamed: 0" in d.columns:
        d = d.rename(columns={"Unnamed: 0":"date"})
    else:
        raise ValueError("未找到日期列（date 或 Unnamed: 0）")
    d["date"] = pd.to_datetime(d["date"])
    for c in d.columns:
        if c != "date":
            d[c] = pd.to_numeric(d[c], errors="coerce")
    return d.set_index("date").sort_index()

def yoy_log(s, k):  # 对数同比（规模类）
    return np.log(s).diff(k)

def yoy_diff(s, k): # 同比差分（占比/比率）
    return s.diff(k)

def one_corr(a, b, method="pearson"):
    # 对齐并计算 r 与 p 值
    idx = a.dropna().index.intersection(b.dropna().index)
    a1, b1 = a.loc[idx], b.loc[idx]
    if len(a1) < 8:
        return np.nan, np.nan, len(a1)
    if method == "pearson":
        r, p = pearsonr(a1, b1)
    else:
        r, p = spearmanr(a1, b1)
    return float(r), float(p), len(a1)

# ===== 读数据 & 同频对齐 =====
emp = read_employment_xlsx(PATH_EMP)
y = emp["employment_rate"].resample(FREQ).mean().rename("y_pred_freq").dropna()

X = read_vars_table(PATH_X)  # 自变量表，已是季度或月度
X = X.resample(FREQ).last()  # 若原本就是同频，这行等效

# 只保留与 y 同期
joint = X.join(y, how="inner")

# ===== 变换：先变换后相关 =====
k = 4 if FREQ == "Q" else 12   # 季度同比滞后4期；月度同比滞后12期

def transform_series(s, name):
    if TRANSFORM == "level":
        return s
    if TRANSFORM == "diff":
        return s.diff(1)
    # 'yoy'
    if "占比" in name:   # 占比做同比差分
        return yoy_diff(s, k)
    else:                # 规模类做对数同比
        return yoy_log(s, k)

y_t = (y if TRANSFORM=="level" else (y.diff(1) if TRANSFORM=="diff" else yoy_diff(y, k))).rename("target")

results = []
for col in joint.columns:
    if col == "y_pred_freq":  # 就业率列已单独处理
        continue
    x_raw = joint[col]
    x_t   = transform_series(x_raw, col)
    r_p, p_p, n_p = one_corr(x_t, y_t, "pearson")
    r_s, p_s, n_s = one_corr(x_t, y_t, "spearman")
    results.append({
        "variable": col,
        "pearson": r_p, "pearson_p": p_p,
        "spearman": r_s, "spearman_p": p_s,
        "n_obs": n_p
    })

tab = pd.DataFrame(results).sort_values("pearson", ascending=False)

# FDR（对 Pearson p 校正）
mask = tab["pearson_p"].notna()
if mask.any():
    _, qvals, *_ = multipletests(tab.loc[mask,"pearson_p"], method="fdr_bh")
    tab.loc[mask,"pearson_q"] = qvals

tab.to_csv(OUT_CSV, index=False)
print(f"已保存：{OUT_CSV}")
print(tab)
